home *** CD-ROM | disk | FTP | other *** search
/ SGI Hot Mix 17 / Hot Mix 17.iso / HM17_SGI / research / lib / igamma_pdf.pro < prev    next >
Text File  |  1997-07-08  |  2KB  |  85 lines
  1. ;$Id: igamma_pdf.pro,v 1.4 1997/01/15 03:11:50 ali Exp $
  2. ;
  3. ; Copyright (c) 1994-1997, Research Systems, Inc.  All rights reserved.
  4. ;       Unauthorized reproduction prohibited.
  5. ;+
  6. ; NAME:
  7. ;       IGAMMA_PDF
  8. ;
  9. ; PURPOSE:
  10. ;       This function computes the incomplete gamma function using a series
  11. ;       representation. It is called by the probability density functions
  12. ;       in this directory. See the function IGAMMA() in the "math"
  13. ;       subdirectory for the user-callable version of the incomplete gamma 
  14. ;       function.
  15. ;
  16. ; MODIFICATION HISTORY:
  17. ;       Modified by:  Jong Yi, Sept 1992
  18. ;                     Increased iterations in g_series.pro
  19. ;       Modified by:  GGS, RSI, July 1994
  20. ;                     Minor changes to code.
  21. ;-
  22.  
  23. pro g_series, result, x, a
  24.   ;Computes incomplete gamma function using a series representation.
  25.   glog = lngamma(a)
  26.   nsample = long(x/50.) > 1000l
  27.   resarray = 1.0/(findgen(nsample) + a)
  28.   resarray[1:*] = x*resarray[1:*]
  29.   sum = 1.0/a
  30.   for i = 1, nsample-1 do begin
  31.     resarray[0] = resarray[0] * resarray[i]
  32.     sum = sum + resarray[0]
  33.     if (abs(resarray[0]) lt abs(sum)*3.0e-7) then begin
  34.       result = sum * exp(-x + a * alog(x) - glog)
  35.       return
  36.     endif
  37.   endfor
  38.   result = -1
  39.   return
  40. end
  41.  
  42. pro g_fract, result, x, a
  43.   glog = lngamma(a)
  44.   gd = 0.0 & fc = 1.0 & b1 = 1.0
  45.   bo = 0.0 & ao = 1.0
  46.   a1 = x
  47.   for n = 1,100 do begin
  48.     an = float(n)
  49.     ana = an -a
  50.     ao = (a1 +ao * ana) * fc
  51.     bo = (b1 + bo *ana) * fc
  52.     anf = an * fc
  53.     a1 = x * ao + anf * a1
  54.     b1 = x * bo + anf * b1
  55.     if a1 then begin
  56.       fc = 1.0/a1
  57.       g = b1 * fc
  58.       if abs((g-gd)/g) LT 3.0e-7 then begin
  59.         result = exp(-x + a * alog(x) - glog) * g
  60.         return
  61.       endif
  62.       gd = g
  63.     endif
  64.   endfor
  65.   result = -1
  66.   return
  67. end
  68.  
  69. function igamma_pdf, a, x
  70.   if x lt a+1.0 then begin
  71.     g_series, result, x, a
  72.     return, result
  73.   endif else begin
  74.     g_fract, result, x, a
  75.     if result ne -1 then return, (1.0 - result) $
  76.     else return, -1
  77.   endelse
  78. end
  79.  
  80.  
  81.   
  82.      
  83.  
  84.  
  85.